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Abstract 



We construct a continuum model for biological aggregations in which individuals 
experience long-range social attraction and short range dispersal. For the case of one 
spatial dimension, we study the steady states analytically and numerically. There 
exist strongly nonlinear states with compact support and steep edges that corre- 
spond to localized biological aggregations, or clumps. These steady state clumps are 
approached through a dynamic coarsening process. In the limit of large population 
size, the clumps approach a constant density swarm with abrupt edges. We use 
energy arguments to understand the nonlinear selection of clump solutions, and to 
predict the internal density in the large population limit. The energy result holds 
in higher dimensions, as well, and is demonstrated via numerical simulations in two 
dimensions. 
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1 Introduction 



Biological aggregations such as insect swarms, ungulate herds, fish schools, 
and bacterial colonies are w i despread exampl e s of self-organization in nature 



( Parrish and Hamnerl. 1997 : Ben- Jacob et al. . 20001 : Okubo and Levin . 2001 



Camazine et all 1200 it ). These groups often arise as social phenomena, without 
direction from a leader or influence of external stimuli such as food and light 
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sources. Social forces among organis ms include attraction, for group cohesion, 



and r epulsion, for collision avoidance (JBreder, 1954; Mo gilner and Edelstein-Keshet 



1999). The resulting aggregati ons can confer benefits such a s prote ction and 



mate choice to their members ( Parrish and Edelstein-Keshet . 1999). 



Mathematical models of social aggregations can be classified into Lagrangian 

and Eulerian types. The Lagrangian approach treats eac h organism as a parti- 

cle obeying a nonlinear d i fference or differential equation (ISakail 1973l:ISuzuki and Sakai , 



1973 [|Okubo et all 



1977 



Vicsek et al.l[l995tlLevine et al. 



2001 



Schweitzer et al 



Erdmann et al.l2002t Parrish et alT 20031 : lAldana and Huepd . 



2001: Couzinetal 



2002: 



2003; lErdmann and EbelinsJ . 120031 : iMorilner et all . [2003). Alternatively, the 



diffusion equation for a continuum p 


3pulation densitv field dKawasakil, |l978l 




Okul 


xl 1980; 


Mimura and Yamaeuti. 


19821: Passo and Demottonil Il984 Hkeda 


1984 


lAltl 11985: Ikcda. 198,4 ISatsuma and Minimal. 198,4 Ikeda and Nasrai 


1987 


iHosono and Mimural.ll989:lGrunbaum and Okubolll994: Edelstein-Keshet et al. 


1998 


iToner and Tul. 1998: 


Flierl et all 1999t Mogilner and Edelstein-Keshetl 




1999 


Topaz and Bertozzil. 20041). A varietv of methods can be used to conned 



the two formulations. The usual method involves a Fokker-Planck approxima- 
tion which relates the distribution of jum p distances made by ind ividuals to 



terms in the advection-diffusion equation ([Okubo and Levinl . 120011 ). Since the 



social communications between organisms often take place a t large dis t ances 



via sight, s ound, or smell, models may be nonlocal in space dKawasaki 
Altl. 11984 llkedal. Il98«4 ISatsuma and Mimural. Il984 Iked a and Nagai l 



1978; 



1987 



1999; 



Hosono and Mimura. Il989: iGrunbaum and OkuboT 19941: Flierl et al- 
Mogilner and Edelstein-Keshet . ll999l : IOkubo et al. . 2001 : Topaz and Bertozzi 
20041) 



Biolog ical aggregations can form distinct groups with sharp edges as they 



move ([Parrish and Hamnerl . Il997t iParrish and Edelstein-Keshetl . 11999). We 



refer to this phenomenon as clumping. While clumping has been observed 



i n two-dimensional numer ical numerical ([Levine et all 120011 ) and analytical 



( Topaz an d Ber tozzil. 2004 ) models, recent mathe matical analyses of swarming 
behaviors ((Mogilner and Edelstein-Keshet , 19991 ) were unable to find the same 
kind of clumping from biologically reasonable one- dimensional models. As we 
will describe below, earlier one-dimensional models which support clumping 
behavior include assumptions about biological interactions which are unlikely 
to be met in nature. 



Our goal in this paper is to investigate how clumping can arise in a sim- 
ple, realistic nonlocal model for biological aggregation, both in one spatial 
dimension and in higher dimensions. Our Eulerian model is an integrodiffer- 
ential conservation law with two movement terms. One describes nonlinear 
degenerate diffusion arising from anti-crowding behavior, and the other de- 
scribes attractive nonlocal social interactions. In contrast to previous studies 
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( Mimura and Yamagutil. 19821: Nagai and Mimura . 19831 : Edelstein-Keshet et al 
19981 Mogilner and Edelstein-Keshet . Il999j ). we do not seek traveling solu- 
tions, which would correspond to cohesive group movement. Rather, we are 
concerned with stationary solutions. Our formulation relates to several earlier 
models that have appeared in the mathematical literature. We discuss these 
now, and then state our main results and outline the remainder of this paper. 



Many nonlocal continuum models for biological aggregation may be cast in 
the form 



dp = d_ /dp 
dt dx \ dx 



VaP 



Here x is the one-dimensional space coordinate, t is time, p(x,t) is the local 
population density, D measures diffusion, and v a is the nonlocal, density- 
dependent speed. 

The earliest models take D constant, and v a = K * u. The asterisk denotes 
spatial convolution and K is an od d spatial weighting function which drops 
off w ith distance and has finite mass ( Kawasaki , 19781 Griinbaum and Okubol . 
1994). K models attractive social forces, which can give rise to spatial instabil- 
ities leading to a unique steady state with a non -uniform spatial distribution 
of the population (|Griinbaum and Okubo, Il994l ). However, these patterns do 
not include clumps (aggregations with compact support). 



Extensions by Mogilner and Edelstein-Keshed (|l999j ) consider group drift by 
including a local, density-dependent velocity as well as by including an even 
component in K, giving rise to traveling swarms. These swarms are not stable 
over long periods of time, ha ving a tendency to break down by losin g individu- 



19991) . However, 



als at the rear of the swarm ((Mogilner and Edelstein-Keshetl . 
the authors note (via analysis and numerical experiments) that nonlinear dif- 
fusion has a tendency to stabilize the swarms. 



On the other hand, earlier e x tensions of the work of Kawasa ril dl978l) by 
Mimura and Yamagutil (Il982h iNagai and Mimural (|l983h . lAltl (|l985h . Ilkeda 
(1985) and Ike da and Naga I (|l987l ) include density dependent (rather than 
co nsta n t) diff usion, still coupled to long-range advective attraction. The model 
of lAltl (|l985h is inspired by chemotactic locomotion, and includes nonlin- 
ear diffusion which is degene r ate for finite p. In contra s t, in the wo rk of 
Mimur a' and Yamagutil (1982). INagai and Mimural (|l983h . Ikeda (1985) and 
Ikeda and Nagail (J1987), the diffusion in (1) takes the form D = ppP~ l . The 
parameter p determines the degree of nonlinearity of the diffusion, and may 
be freely varied. In all of these works, the aggregative weighting function K 
does not meet our previous assumptions of having a finite mass and decaying 
with distanc e. Further, it has support i which is eit her a tunable parameter 
Ikeda (1198,51): Ilkeda and Nagail (|l987l ) or is in finite dMimura and Yamaguti l 
19821 : INagai and Mimural . Il983 l). The model of iHosono and Mimural (Jl989j ) i 



is 
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similar, but adds reaction-type terms modeling logistic growth and preda- 
tion. Depending on y and I, the nonlinear diffusion models hay e stationary 
(llkedal. f l985 ; Ike da and Naea I Il987t lH osono and Mimural . and travel- 



kcua, ryot); ikeaa and iNagai, ryo/; riosono ana Mimura jiygyi i ana travei- 
mg ((Mimura and Yamaguti , Il982l ; Nagai and Mimural . 1 198 31 ) clump solutions. 
Although this class of models exhibits the desirable behavior of clump forma- 
tion, its limitation arises from the biologically unrealistic assumption of strong 
attractive interactions between individuals over arbitrarily large distances. 



In this paper we modify the classic model of Kawasaki ( 1978| ) to include 



density-dependent diffusion. We show that this modification is sufficient to 
give rise to clumped solutions with very sharp edges. These solutions can be 
understood using classical applied mathematical methods of weakly nonlin- 
ear analysis, phase plane analysis, and energy methods. Furthermore, using 
asymptotic and scaling arguments, we show how, in the large population size 
limit, the clumps have constant internal population density. The preferred den- 
sity is predicted by analyzing the ene rgy. The constant internal density p rop- 
erty is typical of biological groups (see Parrish and Edelstein-Keshet , 19991 and 



the extensive discussion in Mogilner et all l2003| ) but has not been found in 



other realistic continuum models. Our model has this property not only in one 
spatial dimension, but in higher dimensions as well. 

The remainder of this paper is organized as follows. We formulate our model 
in Section 2. In Section 3, we outline its basic features, including its con- 
servation properties, linear stability, and energy. In Section 4, we pose the 
model on a periodic one-dimensional domain and perform an in-depth study 
for a particular choice of the aggregative weighting function K. We analyze 
instabilities of trivial population profiles using weakly nonlinear methods. We 
then use phase plane methods to deduce the qualitative structure of swarm 
patterns. We extend the analysis of the qualitative structure through the nu- 
merical calculation of nonlinear steady-state solutions, and the construction of 
numerical bifurcation diagrams. Interestingly, for fixed population size, the al- 
lowed clumps solutions are not unique. We use energy methods discover which 
stationary solutions are preferred, allowing us to deduce the pattern selection 
mechanisms that drive the dynamical system. We also investigate the effects 
of domain size and population size, and perform numerical simulations which 
display the model's coarsening behavior. Section 5 contains a brief extension 
in which we perform numerical simulations in two spatial dimensions to show 
clump formation. Finally, we summarize and conclude in Section 6. 



2 Mathematical model 



Consider a population density p(x, t) > which moves with velocity v(x, t), 
x, v G M™, t > 0. We assume that birth, death, immigration, and emigration 
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of organisms are negligible on the time scale of interest. Then p satisfies the 
standard conservation equation 



pt + V • (vp) = 0. 



(2) 



The members of the population move towards and away from each other fol- 
lowing basic biological principles of aggregation and dispersal, i.e. 



v = v a + v d . 



(3) 



We now develop a kinemat i c mod el ( Edelstein-Keshet et al. , 19981 : Mogilner and Edelstein-Keshet , 
19991 Topaz and Bertozzi . 2004 ) which assumes that the attractive and dis- 
persive velocities v a and depend only on properties of p at the current 
time. 

Though the particulars of aggregation vary from species to species, the sensing 
mechanism responsible (e.g. sight or smell) typically has a characteristic range, 
which we call £, and degrades over distance. For simplicity, we assume that 
the sensing is spatially isotropic, and that organisms sense an averaged nearby 
population. We then associate with an individual in the population at position 
x the sensing function 



s x 



/ K(x - y)p(y) dy = K * p. 



(4) 



The * operator denotes convolution. The kernel K incorporates the sensing 
range and degradation for the particular species under consideration. We as- 
sume, without loss of generality, that K is normalized such that J m „ -ft'(x) c?x = 
£ n . Individuals aggregate by climbing gradients of the sensing function s(x). 
We denote the species-specific characteristic attractive movement speed by V. 
The attractive velocity, then, is 



V£ , 

—V(K*p). 

a 



(5) 



Here, a is a characteristic density. Dimensional arguments dictate that the 
factor of £/at must appear, so that v a has the correct units of velocity. Since 
density has units of number of organisms per unit space, we form the charac- 
teristic density 

a = —. (6) 

That is to say, we take the characteristic density to be the one for which the 
spacing between organisms is the characteristic sensing range. 

Dispersal is assumed to arise a s an anti-crowding mechanism, and operates 



over a much shorter length scale (|Brederl . ll954l : lMogilner and Edelstein-Keshet 



1999). Correspondingly, we take it to be spatially local and in the opposite di- 
rection of population gradients. Moreover, anti-crowding motion is assumed to 
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decrease as the population thins. The simplest model describing these effects 
is linear in both p and Vp, i.e. 



Vri „ 
— ^pVp 

cr 



(7) 



Here, r quantifies the ratio of typical aggregative to repulsive velocities. The 
factor of I /a 2 arises from dimensional arguments, similar to those above for 
the aggregative velocity. 

Combining (2), (3), (5), (6) and (7), we obtain an integrodifferential equation 
for p, namely 



Pt 



V ■ (yt +1 pK * Vp - Vr£ 2n+1 p 2 Vp) = 



where we have used the commutativity of convolution and differentiation. This 
model i s related to the one-dimens io nal models propose d bvlMimura and Yamaguti 
(fmsi. iNagai and Minimal ifmafl. Ifkedal ifmsl and Ifkeda and Nagail (jl987h 
(those models may be put in the same form as ours by using commutativity of 
differentiation and convolution). In those models, the dispersal rate is a (tun- 
able) power of the density, and the sen sing kernels if grow linearly in spac e 



with a possible cutoff at a finite rang e (llkeda L 1985 L Ikeda and Nagail . 
Our model is also similar to that in Mogilner and Edelstein-Keshetl ( 



1987). 



1999). 



However, there, the aggregation and repulsion effects depend directly on the 
same sensing function. Our goal in this paper is to present a model that is 
valid in higher dimensions and incorporates biologically reasonable assump- 
tions about aggregation and dispersal, in particular that the former occurs on 
a longer scale than the latter. 



We rescale the variables as x = (l/£)x, t = (V/£)t, p = (l/a)p and let 
K{x) = K(x). Substituting into (8) and dropping the tildes for convenience, 
we arrive at the dimensionless governing equation 



p t + V • {pK * Vp - rp 2 Vp) = 0. 



(9) 



Due to our rescaling, the characteristic distance of attraction in K is now one, 
and J Rn if(x) <ix = 1. 



As previously discussed, K should reflect tha t biological sensory mechanisms 



have limited spatia l exten t . Common choices ([Mogilner and Edelstein-Keshet 



1999 : ; Okubo et al. . 2001 ) include the decaying exponential and the character- 
istic function: 
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'1/2, 


n = 1 




i^(x) = a n e-W 


a n = < 


1/(2*), 


n = 2 








.l/(4vr 2 ), 


n = 3 








'[-1,1], 




n = 1 




S n = < 


[-1,1] x 


[-1,1], 


n = 2 






.[-1,1] x 


[-1,1] x [-1,1], 


n = 3 



(10a) 



(10b) 



The exponential function arises from assuming a constant rate of transmis- 
sion failure of sensory data per unit distance (i.e. the constant hazard function 
assumption). More complex behavioral models for K% also include differential 
weighting given to stimuli received from different distances. The characteristic 
function K^, on the other hand, arises under the assumption of an identical 
response to all individuals within a fixed distance. 

In order to fully define the model, we must specify a domain D along with 
suitable boundary conditions for (9). We consider three choices. The first is 
no flux boundary conditions, which implies that the members of population 
cannot enter or leave the domain. The second choice is periodic boundary 
conditions, which are mathematically convenient. For this case, the kernel K 
must be modified to be periodic in order for the model to make sense. The final 
choice is to solve the problem in free space, i.e. D = M. n , which corresponds 
to a biological population with no nearby physical barriers. In Section 4 we 
present a case study in one spatial dimension with periodic boundary condi- 
tions. In that section, we discuss how this choice also includes the case of no 
flux boundary conditions, and how it may be used to approximate the free 
space problem. 



3 Basic model characteristics 



The n-dimensional model (9) has a number of useful general properties which 
we describe in this section. In particular, we discuss conservation properties of 
the model and its linear stability. We also show the existence of a Lyapunov 
functional (or energy) which is dissipated under the dynamics of the govern- 
ing equation (9). This is, to our knowledge, the first time such an energy has 
been introduced. We use the energy later to understand the selection of non- 
linear swarm states, and to obtain a quantitative prediction of the equilibrium 
population density for large swarms. 

Conservation of moments. We first discuss the model's conservation prop- 
erties (temporarily restricting attention to the free space problem, i.e. D = 
W 1 ). By construction (see Section 2) (9) conserves the zeroth moment of p, 
namely the total population size or mass M = J D p dx.. Equation (9) also 
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conserves the first moment, or scaled center of mass, M^ l > = f D xp dhc. For 
convenience, we define the inner product (a, b) = J D a • b dx, where the dot 
product must be interpreted appropriately in the case that a and/or b are 
scalars. To show that is conserved, we consider 



M, (1) = (p t ,x) (11a) 

= (-V ■[pK*Vp-rp 2 Vp},x) (lib) 

= (pK * Vp - rp 2 Vp, 1) (11c) 

= (pK * Vp, 1) - (rp 2 Vp, 1) (lid) 

= (pK * Vp, 1) - (rVp 3 /3, 1) (He) 

= (fT*Vp,p)-0 (llf) 

where we have assumed that there is no population contained at infinity. Equa- 
tion (lib) comes from substituting (9), and (11c) and (llf) come from inte- 
gration by parts. We then note that (-£T*Vp, p) = (V[if*p],p) = — (K*p, Vp) 
by commutativity of differentiation and convolution, and integration by parts. 
Also, we have that (K * Vp, p) = (Vp, K * p) = (K * p, Vp) since the convolu- 
tion may be moved across the inner product. Since {K*p, Vp) = —(K*p, Vp), 
= and hence the center of mass is conserved. It follows directly that (9) 
does not support any solution, including unidirectional traveling solutions, in 
which the population undergoes a net translation. Such solutions could exist 
if we included additional advective terms, for instance those due to gradients 
in food sources. We omit such advective terms since we are primarily inter- 
ested in studying the balance of aggregative and dispersive effects. Equation 
(9) does not conserve higher moments. For instance, the spread of the popu- 
lation is free to change under the dynamics of (9), and indeed, will do so as 
the populations clumps (see Section 4.4). 

Linear stability. Equation (9) admits steady states with constant density 
Po > 0. To analyze their stability, we let p = po + pe tqx+<J ' where q is a per- 
turbation wave vector and a is the linear growth rate. Note that the admissi- 
ble q depend on the chosen boundary conditions. For the free space problem, 
q G R n . For the case of periodic boundary conditions, i.e. when D is n-torus of 
size L, q G Z n -2ir/L. For the case of no-flux boundary conditions, q G Z n -7r/L. 
Substituting into (9) and retaining only the terms linear in p, we obtain the 
dispersion relation 

a(q)=p q 2 (K(q)-rp ). (12) 

where q — |q| and K denotes the Fourier transf orm of K. Note that o~(0) = 
for all po due to the conservative structure of (9) ([Cross and Hohenberd . ll993h . 
We determine the neutral stability curve q(po) by setting a = in (12). Figure 
1 shows examples for n = 1 spatial dimension in free space. We consider 
both choices of interaction kernel from (10a), for which K\ = 1/(1 + q 2 ) and 
Kf = (l/g)sin(g). For the latter case, Kf is oscillatory, so there may be 
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Fig. 1. (a) Neutral stability curve for (9) calculated from the dispersion relation a(q) 
in (12) with K = Kf from (10a). The region below the curve corresponds to linear 
instability, (b) Like (a), but with K = Kf from (10a). Since Kf is discontinuous, 
a(q) is oscillatory and the neutral curve has multiple branches. The regions inside 
the curves correspond to linear instability. For both (a) and (b), r = 1 in (9). 

multiple bands of unstable q which are manifest as the multi-branch curve in 
Figure lb. The former case will be considered in Section 4. 



Energy. Finally, we note that (9) possesses a Lyapunov function or energy 

E(p)= I r -p 3 - P K*pdx (13) 
jd 3 

where the first term arises from avoidance and the second from aggregation. 
This energy is dissipated under the dynamics of (9). To see this, we note 



Ek = ( Pt , rp 2 ) - (p t , K*p)-(p,K* Pt ) (14a) 

= {p t , rp 2 ) - (p t , K*p)- (p t , K * p) (14b) 

= ( Pt ,rp 2 )-2(p tl K*p) (14c) 

= {pt, rp 2 - 2K * p) (14d) 

= (-V • [pK * Vp - rp 2 Vp},rp 2 - 2K * p) (14e) 

= (pK * Vp - rp 2 Vp, 2rpVp -2K* Vp) (14f) 

= -2 I p(\K *Vp\-rp\Vp\) 2 rfx < (14g) 

JD 

where (14b) follows from moving the convolution across the inner product, 
(14e) follows from substituting (9), (14f) follows from integration by parts, and 
the inequality in (14g) follows since p > because it is a density. Although 
Lyapuno v function s appear in earlier analysis of one-dimensional swarming 
models (Alt, 1985| ). the particular energy (13) has not been discussed be- 



fore. This non-convex functional is composed of a positive avoidance term 
J D hp 3 g?x and a negative aggregation term — J D pK * p <ix which have differ- 
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ent nonlinear dependence on p and different length scales. This is the hallmark 
of many pattern- forming systems in nature. Non-biological examples include 
the Cahn-Hi lliard equation, which models th e spinodal decomposit ion of bi- 
nary alloys f Cahnl. 11968c lEilbeck et all 1989: Bates and Fife, 1990), droplet 



form a tion in dewetting fluid film s ((Bertozzi et all l20o"l[ " lOron and Bankofi 



2001; Glasner and Witelskil 



2003|) . and self-aggregation of finite-sized parti- 
2005). These models are well-known to exhibit 



cles ((Holm and Putkaradze 
coarsening dynamics, in which small localized clumps form and merge into 
larger clumps over time. We will demonstrate that (9) exhibits similar behav- 
ior. 



4 Case study in one spatial dimension 



We now consider in detail (9) posed in one spatial dimension, i.e. n = 1. Our 
primary goal is to determine the possible steady states for different popula- 
tion sizes, and thus we conduct a bifurcation study with the population size 
M as the bifurcation parameter. We derive analytical expressions for weakly 
nonlinear solutions, and use phase plane, energy, and numerical methods to 
investigate the more nonlinear clumped solutions. We take K = Kf, the expo- 
nential kernel in (10a), because of the biologically reasonable constant hazard 
function assumption that leads to it. For mathematical convenience, we begin 
by considering a periodic domain, which requires that K\ be modified so that 
it is likewise periodic. A simple way to do this is to note that in free space, 
the operator K{* is the inverse of the operator I — d 2 x . Thus, a natural modi- 
fication is to choose a periodic kernel which is precisely the inverse of / — d 2 
on a periodic box. This way, the free-space and periodic kernels have Fourier 
transforms which coincide (or more precisely, the continuous transform of the 
free-space kernel interpolates the discrete transform of the periodic one). In 
Section 4.5, in order to approximate (9) in free space, we take the limit of 
large domain lengths L. Furthermore, we note that to study the case of no- 
flux boundary conditions, one may simply construct an even extension of the 
problem, which then has periodic boundary conditions. 



4-1 Weakly nonlinear analysis 



We pose (9) on a periodic box D = [0, L] of length L = and consider the 
instability of the constant density steady state as the total population size M 
is decreased through the critical value 

M c = L Pc = = *M (15) 
q qr 
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where p c is the critical constant density for the given box size, as determined 
by the dispersion relation (12). We seek weakly nonlinear solutions of (9) via 
a multiple-scales perturbation expansion, letting 



d t ^d t + e 2 d T (16a) 
M = M c + e 2 m (16b) 
p(x,t) = p + ep 1 (x,T) + e 2 p 2 (x,T) + e 3 p 3 (x,T) . . . . (16c) 

Here, e C 1 is a small bookkeeping parameter, T is a slow time scale, and 
m is the deviation from the critical population size. The pi do not evolve 
on the fast time scale t, and so the t dependence is henceforth eliminated 
from our calculation. The perturbation expansion allows us to compute a 
solution to the nonlinear homogeneous problem (9) in terms of the expansion 
functions po(x,T), pi(x,T) and so on, where these functions satisfy linear 
nonhomogeneous equations. Also, 

Po = M/L (17) 

by definition, and 

Pl = Z {T) e iqx +c.c. (18) 

where c.c. stands for complex conjugate and z is the slowly varying amplitude 
of the solution, for which we will derive an amplitude equation. 

At 0(e), (9) becomes the linear problem 

p c K * p ljXX - rp 2 c p hxx = (19) 

which simply recovers the dispersion relation (12). At C(e 2 ), we obtain 



p c K * p 2 , xx - rp 2 c p 2 ,xx = ~pi,xK * p x , x - p\K * Phxx (20) 

+ 2rp c (p lja; ) 2 + 2rp c p 1 p hxx 

which has the particular solution 

p 2 = ^ z 2 (T) e 2iqx +c.c. . (21) 

' 2K{q)-2K{2q) V 1 V ' 

At C(e 3 ), we obtain 



PcK * p^ xx - rp\p z , xx = - p hT - p hx K * p 2;X - p 2 , x K * p hx (22) 

- piK * p 2)XX - (m/L)K * p l)XX 

- p 2 K * p ljXX + Arp c p 1:X p 2:X + 2rp 1 (p 1:X ) 2 
+ 2rp c p 1 p 2 , xx + rpjp 1>xx + 2rp c p 2 pi, xx 

+ 2r(m/L) Pc p hxx . 



11 



The right hand side of (22) contains terms with spatial dependence e tgx which 
lie in the solution space of the left hand side. To guarantee that p is periodic, 
we must eliminate these secular terms. Enforcing this solvability condition 
leads to the amplitude equation 

^- = \ z + a\z\ 2 z, X = -mK(q)q 2 /L, a= _ r ^ K ^ . ( 23 ) 

dT 11 ' Wy 1 ' 2K(q) - 2K(2q) 1 ' 

This result is valid for any K satisfying the assumptions made in Section 2. 



Equation (23) describes a pitchfork bifurcation ( Crawford! . 1991 ) from the con- 



stant density state to a spatially-periodic pattern with steady-state amplitude 
determined by \z\ 2 = —Xm/a. For K = K\ from (10a), a > and so the 
patterned state bifurcates subcritically, and hence unstably. In fact, this is the 
case for any K and q satisfying K(q) > K{2q). 

In Section 4.3 we will compare the result in (23) with numerically obtained 
results. The quantity we will study is the amplitude A = max^p — mini) p. 
Keeping the first two terms in our perturbation expansion, i.e. p = po + epi, 
and substituting (23), we arrive at the weakly nonlinear prediction 



A = -J^(Mr- K{q)L) (K{q) - K(2q)). (24) 

It follows that the amplitude of the periodic pattern is inversely proportional 
to r, or directly proportional to the ratio of typical attractive speed to typical 
repulsive speed. 



4-2 Phase plane analysis for steady states 



We set the flux in (9) equal to and divide by p (discarding the trivial solution) 
to obtain the steady state problem K\ * p x — rpp x = 0. Inverting K\ yields 
the local steady state equation 

Px = rpp x - r(pp x ) xx (25) 

which may be integrated once to obtain p = rp 2 /2 — r(pp x ) x + C/2, where C/2 
is the constant of integration. This equation has an exact solution in terms of 
elliptic integrals, but it is difficult to analyze. Instead, we study the equivalent 
phase-plane problem 

P = <t>, = 



1 p <p 2 C 

r 2 p 2rp 
rp 4 ^ r/ro 2 r ' - a 



4 



(26) 
(27) 
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p=0 p=0 

Fig. 2. (a) - (d): phase portraits for (26) corresponding to different values of C for 
the cases cases I, III, V, VII described in detail in the text. The phase plane problem 
describes the steady states of (9). The stable and unstable manifolds of the fixed 
point F + are indicated by the broken line, and the critical touchdown solution T is 
indicated by the dotted line. 

where = p x , the dot represents d/dx, and £ is the "energy" (different from 
the energy functional in equation 13) which is constant on trajectories of the 
steady state problem (26). We classify the solutions as the parameter C is 
varied. Cases I, III, V, and VII below are depicted in Figure 2. 

Case I: C > 1/r. There are no finite attractors. All trajectories are unbounded. 

Case II: C = 1/r. A fixed point is born at (p, 0) = (1/r, 0). 

Case III: 8/(9r) < C < 1/r. There are two fixed points correspond- 
ing to constant density populations, located at (p, 0) = (p ± ,0), p ± = (1 ± 
Vl — rC)/r. F + is a saddle whose stable/unstable manifold W s ' u forms a ho- 
moclinic loop which is bounded away from the axis. F~ is a nonlinear center, 
and lies inside W s ' u . Between F~ and W s ' u lies a continuous family of periodic 
orbits which correspond to smooth, patterned populations. 

Case IV: C = 8/(9r). W s ' u collides with the origin. 
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Case V: < C < 8/(9r). The homoclinic loop is destroyed, i.e. W s,u no longer 
self-intersects. Periodic orbits are destroyed as they collide with the origin. 
There exists a critical orbit corresponding to a population with peak density 
Pmax — (4 — V16 — 18rC)/3. This touchdown solution just touches down to 
p = and does so with finite slope <p T = ±^/C/(2r). The region bounded by 
the touchdown trajectory, the 0-axis, and W s,u contains a continuous family 
of trajectories which approach (p, <p) = (0, ±oo). For p <C 1 and |0| ^> 1, (26) 
is approximately pp + (p) 2 = 0, which has positive nontrivial solutions of the 
form yax + b. Thus, solutions in S reach p = for finite x and correspond to 
compactly supported swarm-like populations which touch down with infinite 
slope. 

Case VI: C = 0. F~ collides with the origin. 

Case VII: C < 0. The periodic orbits are destroyed. The only bounded solu- 
tions are the compactly supported ones, lying between W s ' u and the 0-axis. 

In summary, for 8/(9r) < C < 1/r, only trivial (constant density) and periodic 
solutions are possible. For < C < 8/(9r), trivial, periodic, and clumped 
solutions are possible. For C < 0, only trivial and clumped solutions are 
possible. Note that the solutions of (26) actually form a two parameter family. 
One parameter is the mathematical constant of integration C, which selects 
a particular phase portrait. The other parameter may be taken to be the 
"initial condition" in that phase portrait, which selects a particular trajectory. 
More biologically relevant parameters include the total population size M, the 
peak density p max , and the period or support of the steady solution. In fact, 
fixing any two of these quantities determines the third. Unfortunately, the 
relationship between these biological quantities and the mathematical ones we 
have used for our phase plane study cannot be determined by the qualitative 
analysis we have carried out in this section. In order to understand the role 
of the biologically relevant quantities, we perform a numerical investigation in 
the next section. 



4-3 Numerical solution and bifurcation diagram 

We fix L and construct the bifurcation diagram describing the steady state 
solutions of (9). An example for L = 2tt with r = 1 is given in Figure 3a, which 
shows the solution amplitude A = max D p — min D p as the population size M 
is varied. Solution profiles corresponding to different values of M along the 
bifurcation curve are sown in Figure 3b. The bifurcation diagram and solution 
profiles are obtained in the following way. 

For the periodic and touchdown solutions, we compute trajectories of (26) 
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Fig. 3. (a) Bifurcation diagram for (9) with K = K e from (10a) depicting solution 
amplitude A versus the total population size M. The size of the periodic box is 
L = 2tt, and we take r = 1. The different solution types are stable constant den- 
sity (horizontal sold line), unstable constant density (dashed line), periodic (open 
circles), touchdown periodic (x), and compactly supported (dots). The solid curve 
is the bifurcation curve predicted by (23), the result of a weakly nonlinear analysis, 
(b) Steady state solutions for different values of M corresponding to different points 
along the bifurcation curve in (a). 



numerically using the Runge-Kutta (4,5) method, locate the trajectory with 
period L, and measure M (using the trapezoidal rule) and A. This process is 
performed different values of C. We perform a similar procedure for the clump 
solutions. However, for fixed C, any clump solution with the correct mass and 
with support that fits inside D is allowed. Thus, for fixed M, there is not a 
unique clump, but rather a continuous one-parameter family, of which each 
member has a different length and peak density. There is a selection mechanism 
which prefers one member of this family for each M (in particular, the member 
represented in Figure 3a). We discuss this mechanism in Section 4.5. 

Figure 3a reveals that the subcritical branch of periodic solutions connects 
to the clump solutions at the touchdown solution. The subcritical branch of 
clump solutions turns around, and eventually terminates at the empty pop- 
ulation solution M — 0. From the bifurcation diagram we see that as M is 
decreased through the bifurcation point and the trivial state destabilizes, the 
transition will involve a jump to a large-amplitude swarm. One interesting bi- 
ological prediction which follows is that of hysteretic transitions for organisms 
in a confined environment (in which case the appropriate no-flux boundary 
conditions may be related to the periodic ones we use here; see previous dis- 
cussion). Large, homogeneous populations will not form clumps. If sufficiently 
many organisms are removed, so as to bring the population size below the 
critical one, clumping will occur. Then, if organisms are added back in, the 
population may remain clumped even for population sizes above the critical 
one (so long as not too many are added). 
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Fig. 4. Peak density p m ax versus population size M of clump solutions as computed 
on periodic domains of length L = 2ir (dots) and 4ir (x's), both with r = 1. The 
two data sets coincide for small M, but diverge for larger M when "wraparound" 
effects due to periodicity come into play. The L = 2tt branch terminates for smaller 
M since compactly supported solutions with larger M do not fit in the domain. 
The dotted vertical line indicates the critical population size for L = 2tt, namely 
M c = vr (cf. Figure 3). 

The particular form of the connection to the connection of the periodic branch 
to the branch of clumps is an artifact of having a finite periodic domain. As 
discussed in Section 2, for the periodic problem, the kernel K\ must be made 
periodic, and hence depends on the domain length L. Therefore, steady state 
solutions of identical mass, but on different domains, will not necessarily be 
identical. Figure 4 demonstrates this effect. We plot the maximum density 
Pmax versus population size M of compactly supported solutions as computed 
on periodic domains of length L = 2n and 4tt. The two results agree for 
small M, but diverge for larger M. One intuitive way of understanding this 
phenomenon is to imagine that for clumps with sufficiently large support, the 
left and right sides of the clump may interact strongly due to the nonlocal 
aggregative term in (9) and the periodic boundary conditions. 



4-4 Coarsening 



It is of interest to investigate how the steady state clumps just discussed 
are approached. Numerical simulations reveal that our model displays "coars- 
ening" dynamics. We use a fully implicit numerical scheme. Derivatives are 
calculated using central finite differences, and the convolution operator Kf* 
is computed as (J — c^) -1 , as discussed previously. The nonlinear problem is 
solved via Newton iteration. 

As an example, we perform simulations on a periodic box of length L = 8ir 
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Fig. 5. Snapshots showing the coarsening of an initially randomly-distributed pop- 
ulation with M = 10 on a domain of length L = 8n. Here, r = 1 in (9). 

and choose a random initial condition with total integral M = 10. Figure 5 
shows a series of snapshots of the solution as it is evolved according to (9). 
The irregular initial profile is rapidly smoothed. Then a more slow transition 
occurs, in which the smooth profile breaks up, in this case into three distinct 
clumps of organisms. However, due to the long-range social attraction modeled 
in K, each clump can feel the presence of the others. The two rigthmost clumps 
move slowly towards each other and merge into a larger group. Then, this large 
group and the remaining group are attracted together, though the motion 
occurs over an even slower time scale. Eventually, the two groups merge into 
one large group, and steady state is reached. 



The coarsening dynamics are of special interest to us given the widely held view 
that "social behaviors that on short time scales lead to the formation of social 
groups. . . at the largest time and space scales [have] profound consequences 
for ecosystem dynamic s and for the evolu tion of behavioral, morphological, 



and life history traits" (O kubo et all 120011 ). The vast majority of models for 



splitting and joining of social groups on long time and space scales are stochas- 
tic, and involve quantities such as the probability of a group splitting into k 
smaller groups, which may be extremely difficult to measure experimentally. 
Coarsening dynamics are a natural context in which to study the splitting and 
merging of groups. Moreover, this formulation has the advantage that results 
(such as the coarsening rate) could be explicitly tied to the rules for movement 
(for example, the properties of the interaction kernel K, or the value of the 
speed ratio r). 



4-5 Nonlinear selection and large-domain, large-population asymptotics 



As discussed in Section 4.2, for fixed population size M there is a continuous 
family of clump solutions. However, as stated previously, and as demonstrated 



17 




-0.89 



(b) 



-0.9 



E 



-0.91 



x 



-0.92 



0.7 



0.8 



0.9 



x 



P, 



Fig. 6. (a) Members of the continuous family of compactly supported steady state 
solutions to (9) having total population size M = 2.51. Solutions are obtained 
from numerical integration of the phase-plane problem (26). The thick broken line 
corresponds to the profile having the minimum energy E in (13), which is the one 
observed in simulations of (9). (b) Energies of the solutions in (a), parameterized 
by maximum density p max - The minimum is indicated by an "x". For this example, 
r = 1 in (9). 

by the simulations of Section 4.4, one member of this family is preferred by the 
dynamics of (9), and it is precisely the one which minimizes the energy (13). 
Figure 6 verifies this selection mechanism. Figure 6a shows some members 
of the continuous family, calculated from (26) for M = 2.51 with r = 1 in 
(9). The minimum-energy solution appears as a thick broken line. Figure 6b 
shows the energy E of the solutions, parameterized by their maximum density 
Pmax- To verify the energy argument, we evolve a random initial condition with 
M = 2.51 according to (9) using an implicit, central in space (1, 2)-accurate 
finite difference method. When plotted on Figure 6a, the steady state achieved 
in the simulation is indistinguishable from the energy-minimizing solution in 
6a calculated from (26). 

Though thus far we have considered small periodic domains, it is perhaps 
more biologically relevant to consider the steady states of (9) for the free 
space problem, i.e. with D = 1L We approximate these solutions by taking 
an extremely large domain length L. We calculate the clump solutions for 
different population sizes M with L taken to be large enough for each M so 
that the previously discussed finite-domain effect is negligible (a rule of thumb 
is that for a fixed L, we calculate only those clumps with M < 0.5M C ). In 
Figure 7a, we plot p max as a function of M. Selected density profiles are shown 
in Figure 7b. For small values of M, e.g. M = 2 and M — 7, the clumped 
solutions have small p max and a narrow, rounded shape. For larger population 
sizes, e.g. M = 20, p max is around 1.5, and the shape of the profile is still 
quite rounded. For much larger population sizes, e.g. M = 50 and M = 80, 
the clumps approach a rectangular shape, i.e. they have a nearly constant 
internal density, deviating only very close to the edges of the group which go 
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Fig. 7. Peak density p m ax versus population size M of clump solutions computed 
on periodic domains of very large length L so as to approximate solutions on K (see 
text for details). In the large M limit, density profiles have p ma x = 1-5 as predicted 
by the energy (13). (b) Density profiles for various values of M. Only the center of 
the domain, which contains the entire support of the clumps, is shown. 



steeply to p = 0. Interestingly, as M is made larger and larger, p max stays at 
1.5 (see Figure 7a) and the rectangular profile simply becomes wider. 

The asymptote at p max = 1.5 in Figure 7a may be understood from scaling and 
energy arguments. If we assume that the steady solution has a characteristic 
density p* and a characteristic length scale L*, then the left side of the steady 
state equation (25) scales as p^/L*, the first term on the right hand side 
scales as pi/ L*, and the remaining term scales as pi/ ' L\. In the limit of "long" 
solutions, the balance is between the left hand side and the first term on 
the right hand side, with the remaining term being a much smaller 0(1/ L 2 ) 
correction. Then we have p*/L* ~ pl/L*, and so p* is an 0(1) constant. If 
we assume the population profile approaches a rectangular shape, we may use 
the energy (13) to determine the height p* of the rectangle. For rectangular 
solutions, (13) yields 



E = I T -pl- P*K * p* dx. (28) 
Jd 3 



K has characteristic length scale 1 (see Section 2), and so for sufficiently large 
L*, K resembles a 5-function and convolution with K acts as the identity. 
Then we have 
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(29d) 



(29b) 



(29a) 



(29c) 



where we recall that for a rectangular solution, M = p*L*. As previously dis- 
cussed, for fixed M, the dynamically preferred steady state is the one which 
minimizes E. Minimizing (29d) with respect to p*, we find p* = 3/(2r), in 
agreement with our numerical observations for r — 1. From the inverse de- 
pendence of p* on r, we have the general result that the constant density of 
large rectangular groups is directly proportional to the ratio of characteristic 
attractive to repulsive velocities. Note that all of the above arguments are 
independent of the particular choice of interaction kernel K after rescaling. 



5 Simulations in two spatial dimensions 

The energy argument given above is, in fact, independent of the number of 
dimensions, and so the preferred density for large clumps is always expected 
to be p* = 3/(2r). This is demonstrated by numerical simulations. We conduct 
our simulations on a periodic box. We use a hybrid numerical method with 
adaptive time stepping. The attraction term is treated explicitly and pseu- 
dospectrally. The dispersal term is treated implicitly, with operator splitting. 
For each time iteration, dispersal in the x and y directions are implemented 
successively (the order is alternated). The derivatives are approximated using 
central finite differences, and the resulting nonlinear problem is solved with 
Newton iteration. 

Figure 8 shows a numerical simulation of (9) in a two dimensional box with 
each side L = 40. The population size is M = 600 and the speed ratio r — 1. 
We take N = 300 points along each axis. The initial condition (Figure 8a) con- 
sists of two disjoint discs, with a randomly distributed population in each one. 
The density profile within each disk rapidly smooths out (Figure 8b). Over a 
slower time scale, the two discs merge (Figure 8c). The resulting population 
clump retains a thin interfacial region outside of which p = 0. The interface 
evolves on a very slow time scale, and the internal population density becomes 
nearly constant (Figure 8d). Because the evolution of the group boundary is 
very slow, we terminate the simulation before steady state is reached. Nonethe- 
less, Figure 9 shows the peak density p max as a function of time, and confirms 
that it approaches 1.5, as predicted by the energy arguments. 
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Fig. 8. (a-d) Population density profiles under the dynamics of (9). The initial con- 
dition consists of two disjoint discs, with the population randomly distributed in 
each. The population in each disc rapidly smooths out, and the two discs merge 
over a slower time scale, and form a group with a nearly constant internal popu- 
lation density. The simulations are conducted with population size M = 600 and 
characteristic repulsive to attractive speed ratio r = 1 on a box of size L = 40 with 
N = 300 gridpoints on each axis. 

6 Conclusion 

The purpose of this paper is to demonstrate how population clumping can 
occur in a realistic, minimal model for the movement of organisms. We have 
focused on a kinematic model, i.e. one in which the velocity is a functional of 
the density (as opposed to a dynamic grouping model, which would include 
inertial forces). We have shown that the balance between nonlinear diffusion 
and nonlocal social attraction (which decays with distance) is sufficient for 
the formation of localized groups from initially disperse population profiles. 
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Fig. 9. Peak density p ma x as a function of time for the simulation shown in Figure 8. 
The peak density approaches 1.5, as predicted by energy arguments. 



Our focus h as been on stationary , rath er than traveling, solutions. Following 
the ideas in iTopaz and Bertozzil (|2004l ). group motion may be incorporated 
into the model by incorporating additional social interaction terms of incom- 
pressible form, which have been shown to provide cohesive group motion while 
preserving the properties of constant density and compact support. Alterna- 
tively, group motion may be provided by considering external factors such as 
movement up food gradients. Of course, our model can also be interpreted 
simply as one for the formation of stationary clumps of organisms. 



Our case study in one spatial dimension uncovered three key properties which 
agree with biological observation: 



(1) The groups are truly localized, i.e. they have well defined boundaries 
outside of which the population density is 0. 

(2) The density drops steeply to at the edge, i.e. the density profile has 
infinite slope at the boundary. 

(3) For sufficiently large populations, there is a preferred density (or in the 
parlance of discrete models, a preferred inter-organism spacing) indepen- 
dent of the population size itself, i.e. adding more members to the group 
does not change the "packing" of organisms, but simply increases the 
spatial extent of the group. 



Furthermore, for the case of two spatial dimensions, simulations revealed so- 
lutions which also have compact support, steep edges, and a constant inter- 
nal population density. Finally, we have connected a macroscopic property of 
the population clumps, namely the preferred internal population density, to 
a property of the movement rules, namely the ratio of typical dispersive to 
attractive speeds, by analyzing an energy associated with the aggregation. 
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